Stata：最大似然估计入门教程(MLE)-ml
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
连享会 · 2022 空间计量专题
作者：曹昊煜 (兰州大学)
邮箱：caohy19@lzu.edu.cn
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
目录
1. 最大似然估计与条件最大似然估计
1.1 M 估计量
1.2 最大似然估计
1.3 条件最大似然估计
2. 最大似然估计两种方法
2.1 最大似然估计：极简示例
2.2 最大似然估计：网格搜索
3. Stata 范例
3.1 范例一：估计正态分布的参数
3.2 范例二：线性模型
3.3 范例三：定性反应模型
3.4 范例四：断尾回归
3.5 范例五：截取回归
4. 结束语
5. 相关推文
1. 最大似然估计与条件最大似然估计
1.1 M 估计量
计量经济学中一类重要的估计量称为 “极值估计量”，该估计量通过最大化某个目标函数获得，即：
极值估计量又可以分为 M 估计量和 GMM 估计量。最大似然估计属于 M 估计量的一种，其目标函数的一般形式为样本均值：
1.2 最大似然估计
最大似然估计主要是将目标函数中的 替换为密度函数的对数形式。令 为独立同分布的序列，且概率密度为 ，则数据 的联合密度为：
当使用参数的假象值替代真实参数时，联合密度称为似然函数。真实参数 的 ML 估计是选择相应的 ，使得似然函数最大化。由于对数变换是单调的，最大化似然函数等价于对数似然函数最大化：
因此 的 ML 估计量为 M 估计量：
在一定的正则假定下，ML 估计量为一致渐进正态估计量。
1.3 条件最大似然估计
在一般的研究问题中，数据 可以区分为解释变量和被解释变量，而参数估计的主要目标是考察解释变量 对被解释变量 条件分布的影响。根据条件概率密度公式，数据的联合密度可以表示为：
对等式两边取对数，并使用假想的参数替代真实参数得到对数似然函数：
若条件概率密度和解释变量的概率密度之间不存在函数相关，忽略等式右侧的第二项得到简化的对数自然函数。对应的 M 估计量的形式为：
因此该方法称为条件极大似然估计。在一般的应用中，条件极大似然估计更为普遍，其原因在于我们对联合分布中的参数建模后，被解释变量的分布会依赖于解释变量的变化。
2. 最大似然估计两种方法
2.1 最大似然估计：极简示例
首先使用一个简单的例子来说明 MLE，该示例摘自 B 站 Up 主「小崔说数」的讲解内容。假设城市居民只有两种出行方式，分别是步行 () 或开车 ()。二者发生的概率分别为 和 。随机询问五个人，结果为 ，从而可以得到：
| 1 | 0 | 0 | 1 | 1 | |
|---|---|---|---|---|---|
在该模型中，相当于估计参数 。由于样本是随机抽取的，因此联合概率密度等于所有样本概率之积，定义似然函数：
对等式两边取对数，则有：
求取一阶条件：
估计结果 与 的样本均值相同，因为在二值观测中，期望与取值为 1 的概率相同。通过上面这个简单的例子，可以看出 MLE 的基本步骤：
第一步：确定 的概率分布； 第二步：写出样本中所有观测值的联合分布 ； 第三步：极大化 从而得到 的估计。
2.2 最大似然估计：网格搜索
求解最大似然估计 (极大似然估计) 的结果可以使用一阶条件，或称为似然方程，但在某些情形下，并不存在参数的显式解，因此需要使用数值方法求解参数。一种简单的，便于考察 MLE 整个过程的方式是网格搜索，当真实参数范围已知时，可以通过搜索的方式估计参数。其过程如下：
以某个特定的步长在参数的取值范围内选择不同的参数值； 计算特定参数取值下的似然函数值； 选取最大化似然函数取值对应的参数作为真实参数的估计。
最大似然估计的一个简单例子是估计正态分布的均值和标准差。由正态分布的概率密度函数，平均对数似然函数为：
假定数据生成过程服从均值为 1，标准差为 0.8 的正态分布，则网格搜索的基本思路是同时选择不同的 和 ，并计算对数似然值并比较其大小。过程如下：
. clear
. set obs 10000
. set seed 1010110 // 设定种子值
. gen y = rnormal(1,0.8) in 1/1000 // 生成随机样本 y ~ N(1,0.8)
. gen mu = .
. gen sigma = .
. gen logL = . // 生成用于保存均值、标准差和对数似然值的
. local n = 1
. * 开始网格搜索
. forvalue mu = 0.01(0.01)1{
2. forvalues sigma = 0.01(0.01)1{
3. gen term1 = (y-`mu')^2/(2*`sigma'^2)
4. egen term2 = sum(term1)
5. qui sum term2
6. local last = r(max)
7. local N = 1000
8. replace logL = -0.5*`N'*log(2*_pi) - 0.5*`N'*log(`sigma'^2) - `last' in `n' // 保存对数似然函数值
9. replace mu = `mu' in `n'
10. replace sigma = `sigma' in `n'
11. drop term1 term2
12. local n = `n'+1
13. }
14. }
. * 搜索结果
. gsort - logL
. list mu sigma logL in 1/10
+-------------------------+
| mu sigma logL |
|-------------------------|
1. | .99 .81 -1208.978 |
2. | .99 .82 -1209.109 |
3. | .99 .8 -1209.152 |
4. | .99 .83 -1209.527 |
5. | .99 .79 -1209.652 |
|-------------------------|
6. | .98 .81 -1209.714 |
7. | .98 .82 -1209.827 |
8. | .98 .8 -1209.907 |
9. | .99 .84 -1210.215 |
10. | .98 .83 -1210.228 |
+-------------------------+可以看出，均值和标准差的搜索结果与真实的参数非常接近。当然，在具体的 MLE 中，我们并不知道真实参数的取值范围，因此需要使用牛顿-拉夫森等其他数值计算方法，在以下的例子中 Stata 会自动实现数值计算。
3. Stata 范例
在 Stata 中，可以使用 ml 命令灵活地实现各类最大似然估计。关于 ml 命令的使用和设定方式可以参考帮助文档 help ml 或连享会推文「Stata：数值求解极大值及 MLE 示例」。本文将通过五个简单的范例，介绍如何在 Stata 中实现基本的 MLE 方法。
3.1 范例一：估计正态分布的参数
第一个例子仍然是对正态分布中均值和标准差的估计。此时我们不再使用网格搜索，而是使用其他数值方法，概率密度函数仍然服从正态分布：
在 Stata 中使用 ml 命令，只需要定义特定观测值下的概率密度，而不是整体的对数似然函数，例如观测值 的概率密度是均值为 ，标准差为 的正态分布概率密度函数。
. clear
. set obs 1000
. set seed 1010110
. gen y = rnormal(1,0.8) // 生成随机样本 y ~ N(1,0.8)
. cap program drop mynormal_lf
. program define mynormal_lf // 定义对数似然函数
1. version 16
2. args lnf mu sigma
3. qui replace `lnf' = ln(normalden($ML_y1,`mu',`sigma'))
4. end
. ml model lf mynormal_lf (y = ) (sigma: ) // 模型形式
. ml maximize // 参数求解
Iteration 0: log likelihood = -1247.1045
Iteration 1: log likelihood = -1207.7554
Iteration 2: log likelihood = -1207.5499
Iteration 3: log likelihood = -1207.5496
Iteration 4: log likelihood = -1207.5496
Number of obs = 1,000
Wald chi2(0) = .
Log likelihood = -1207.5496 Prob > chi2 = .
------------------------------------------------------------------------------
y | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
eq1 |
_cons | 1.033 0.026 40.37 0.000 0.983 1.083
-------------+----------------------------------------------------------------
sigma |
_cons | 0.809 0.018 44.72 0.000 0.774 0.845
------------------------------------------------------------------------------关于似然函数的定义需要注意以下几点：
cap program drop：用于删去可能存在的重复程序；progname_lf：是程序命名的一般形式，即自己确定的程序名称加上下划线和lf后缀；version 16：声明 Stata 版本；args：后面紧跟对数似然函数值lnf和参数名称；$ML_y1：代表被解释变量。
估计结果中显示了 MLE 的迭代过程，最大化的对数似然值为 -1207.5496，估计出的均值为 1.033，标准差为 0.809，与设定的数据生成过程基本相同。
3.2 范例二：线性模型
对于均值和标准差的估计属于无条件最大似然估计，因为我们并没有对分布中的参数进行建模，也相当于对截距项回归。接下来考察普通的线性模型：
假设随机干扰项 。从而 ，即对正态分布的均值进行了线性建模。相应的条件对数似然函数为：
使用 Stata 进行线性模型的 ML 估计，并与 regress 命令做对比：
. sysuse auto.dta, clear
. cap program drop mylinear_lf
. program define mylinear_lf
1. version 16
2. args lnf xb sigma
3. qui replace `lnf' = ln(normalden($ML_y1,`xb',`sigma'))
4. end
. ml model lf mylinear_lf (price = foreign trunk weight length) (sigma: )
. ml maximize, nolog
Number of obs = 74
Wald chi2(4) = 90.07
Log likelihood = -666.25274 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
price | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
eq1 |
foreign | 3580.051 624.022 5.74 0.000 2356.990 4803.111
trunk | -10.282 78.824 -0.13 0.896 -164.774 144.210
weight | 5.769 0.934 6.18 0.000 3.938 7.600
length | -89.653 34.534 -2.60 0.009 -157.339 -21.968
_cons | 4672.931 3852.845 1.21 0.225 -2878.506 12224.368
-------------+----------------------------------------------------------------
sigma |
_cons | 1967.417 161.721 12.17 0.000 1650.451 2284.384
------------------------------------------------------------------------------
. reg price foreign trunk weight length
Source | SS df MS Number of obs = 74
-------------+---------------------------------- F(4, 69) = 21.00
Model | 348631330 4 87157832.5 Prob > F = 0.0000
Residual | 286434066 69 4151218.35 R-squared = 0.5490
-------------+---------------------------------- Adj R-squared = 0.5228
Total | 635065396 73 8699525.97 Root MSE = 2037.5
------------------------------------------------------------------------------
price | Coefficient Std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
foreign | 3580.051 646.236 5.54 0.000 2290.845 4869.256
trunk | -10.282 81.630 -0.13 0.900 -173.129 152.565
weight | 5.769 0.967 5.96 0.000 3.839 7.699
length | -89.653 35.763 -2.51 0.015 -160.999 -18.307
_cons | 4672.931 3989.999 1.17 0.246 -3286.899 12632.762
------------------------------------------------------------------------------可见 MLE 与 OLS 的估计结果完全相同，仅仅在标准误上存在细微的差异，如果考察二者的一阶条件，则可以发现二者面临的优化问题是相同的。
3.3 范例三：定性反应模型
接下来考察几个非线性模型的例子。当被解释变量取值为 0-1 时，需要构建定性反应模型，其基本的概率分布为：
对参数 进行建模，由于其代表的是 取 1 的概率，因此其函数形式可以定义为正态分布或 Logistic 分布。当假定 服从正态分布时，该模型称为 Probit 模型：
当假定 服从 Logistic 分布时，该模型称为 Logit 模型：
其中，
在 Stata 中，normal 和 invlogit 分别表示正态分布累积分布函数和拟 Logit 变换，probit 和 logit 命令分别用于估计上述模型。
. * probit 模型
. sysuse auto.dta, clear
. cap program drop myprobit_lf
. program define myprobit_lf
1. version 16
2. args lnf xb
3. qui replace `lnf' = ln(normal(`xb')) if $ML_y1 == 1
4. qui replace `lnf' = ln(1 - normal(`xb')) if $ML_y1 == 0
5. end
. ml model lf myprobit_lf (foreign = price trunk weight length)
. ml maximize, nolog
Number of obs = 74
Wald chi2(4) = 15.11
Log likelihood = -17.72839 Prob > chi2 = 0.0045
------------------------------------------------------------------------------
foreign | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
price | 0.001 0.000 3.32 0.001 0.000 0.001
trunk | 0.010 0.083 0.12 0.901 -0.153 0.174
weight | -0.004 0.002 -2.71 0.007 -0.007 -0.001
length | 0.033 0.045 0.74 0.460 -0.055 0.121
_cons | 1.201 5.070 0.24 0.813 -8.735 11.138
------------------------------------------------------------------------------
. probit foreign price trunk weight length, nolog
Probit regression Number of obs = 74
LR chi2(4) = 54.61
Prob > chi2 = 0.0000
Log likelihood = -17.72839 Pseudo R2 = 0.6063
------------------------------------------------------------------------------
foreign | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
price | 0.001 0.000 3.32 0.001 0.000 0.001
trunk | 0.010 0.083 0.12 0.901 -0.153 0.174
weight | -0.004 0.002 -2.71 0.007 -0.007 -0.001
length | 0.033 0.045 0.74 0.460 -0.055 0.121
_cons | 1.201 5.070 0.24 0.813 -8.735 11.138
------------------------------------------------------------------------------. * logit 模型
. cap program drop mylogit_lf
. program define mylogit_lf
1. version 16
2. args lnf xb
3. qui replace `lnf' = ln(invlogit(`xb')) if $ML_y1 == 1
4. qui replace `lnf' = ln(1 - invlogit(`xb')) if $ML_y1 == 0
5.
. * qui replace `lnf' = ln(exp(`xb')/(1+exp(`xb'))) if $ML_y1 == 1
. * qui replace `lnf' = ln(1 - exp(`xb')/(1+exp(`xb'))) if $ML_y1 == 0
. end
. ml model lf mylogit_lf (foreign = price trunk weight length)
. ml maximize, nolog
Number of obs = 74
Wald chi2(4) = 12.61
Log likelihood = -17.786494 Prob > chi2 = 0.0134
------------------------------------------------------------------------------
foreign | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
price | 0.001 0.000 3.14 0.002 0.000 0.002
trunk | 0.009 0.142 0.06 0.949 -0.269 0.287
weight | -0.007 0.003 -2.55 0.011 -0.013 -0.002
length | 0.050 0.081 0.61 0.539 -0.109 0.209
_cons | 3.388 9.308 0.36 0.716 -14.855 21.631
------------------------------------------------------------------------------
. logit foreign price trunk weight length, nolog
Logistic regression Number of obs = 74
LR chi2(4) = 54.49
Prob > chi2 = 0.0000
Log likelihood = -17.786494 Pseudo R2 = 0.6050
------------------------------------------------------------------------------
foreign | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
price | 0.001 0.000 3.14 0.002 0.000 0.002
trunk | 0.009 0.142 0.06 0.949 -0.269 0.287
weight | -0.007 0.003 -2.55 0.011 -0.013 -0.002
length | 0.050 0.081 0.61 0.539 -0.109 0.209
_cons | 3.388 9.308 0.36 0.716 -14.855 21.631
------------------------------------------------------------------------------从估计结果中，无论是 Probit 还是 Logit 模型，直接编写的最大似然估计结果与官方命令提供的结果基本没有差异。
3.4 范例四：断尾回归
有时候，我们获取的样本并非从总体中随机抽样得到，而是从总体的子样本中获取的，例如规模以上的企业调查。假设随机变量 的分布为 ，则 的断尾分布密度函数为：
假设无断尾 的分布为 ，则在 处断尾的分布为：
从而对数似然函数为：
在 Stata 中，truncreg 命令可以用于估计断尾回归，因此可以使用该命令验证自己编写的断尾回归代码。
. lxhuse laborsub.dta, clear
. keep if whrs > 0 // 假设在 0 处断尾
. cap program drop mytrunc_lf
. program define mytrunc_lf
1. version 16
2. args lnf xb sigma
3. qui replace `lnf' = ln(((1/`sigma')*normalden(($ML_y1-`xb')/`sigma'))/(1-normal((0-`xb')/`sigma')))
4. end
. ml model lf mytrunc_lf (whrs = kl6 k618 wa we) (sigma: )
. ml maximize, nolog
Number of obs = 150
Wald chi2(4) = 10.05
Log likelihood = -1200.9157 Prob > chi2 = 0.0395
------------------------------------------------------------------------------
whrs | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
eq1 |
kl6 | -803.004 321.361 -2.50 0.012 -1432.860 -173.148
k618 | -172.875 88.729 -1.95 0.051 -346.780 1.031
wa | -8.821 14.368 -0.61 0.539 -36.983 19.341
we | 16.529 46.504 0.36 0.722 -74.617 107.674
_cons | 1586.260 912.354 1.74 0.082 -201.922 3374.441
-------------+----------------------------------------------------------------
sigma |
_cons | 983.726 94.443 10.42 0.000 798.621 1168.830
------------------------------------------------------------------------------
. truncreg whrs kl6 k618 wa we, ll(0) nolog
Truncated regression
Limit: Lower = 0 Number of obs = 150
Upper = +inf Wald chi2(4) = 10.05
Log likelihood = -1200.9157 Prob > chi2 = 0.0395
------------------------------------------------------------------------------
whrs | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
kl6 | -803.004 321.361 -2.50 0.012 -1432.861 -173.147
k618 | -172.875 88.729 -1.95 0.051 -346.781 1.031
wa | -8.821 14.368 -0.61 0.539 -36.983 19.341
we | 16.529 46.504 0.36 0.722 -74.617 107.674
_cons | 1586.260 912.355 1.74 0.082 -201.923 3374.442
-------------+----------------------------------------------------------------
/sigma | 983.726 94.443 10.42 0.000 798.621 1168.831
------------------------------------------------------------------------------3.5 范例五：截取回归
截取由 Tobin (1958) 首次引入经济学，与断尾分布不同，此时样本仍然是从总体中随机抽取的，而当被解释变量的取值小于某个特定的值时，全部被归并到该取值上。截取模型又称为归并模型或 Tobit 模型。其形式为潜变量模型：
可以将 的分布视为两个部分， 的部分是正常的观测值，其分布与潜变量 相同，而 的概率与 的概率相同。因此 Tobit 模型的概率密度由两个部分组成：
相应地，条件对数似然函数为：
使用 tobit 命令来验证自己编写的 MLE 程序：
. lxhuse womenwk.dta, clear
. cap program drop mytobit_lf
. program define mytobit_lf
1. version 16
2. args lnf xb sigma
3. qui replace `lnf' = ln((1/`sigma')*normalden(($ML_y1-`xb')/`sigma')) if $ML_y1 != 0
4. qui replace `lnf' = ln(normal((0-`xb')/`sigma')) if $ML_y1 == 0
5. end
. ml model lf mytobit_lf (lwf = age married children education) (sigma: )
. ml maximize, nolog
Number of obs = 2,000
Wald chi2(4) = 472.59
Log likelihood = -3349.9685 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
lwf | Coefficient Std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
eq1 |
age | 0.052 0.006 9.08 0.000 0.041 0.063
married | 0.484 0.104 4.68 0.000 0.281 0.687
children | 0.486 0.032 15.33 0.000 0.424 0.548
education | 0.115 0.015 7.62 0.000 0.085 0.145
_cons | -2.808 0.263 -10.67 0.000 -3.324 -2.292
-------------+----------------------------------------------------------------
sigma |
_cons | 1.873 0.040 46.80 0.000 1.794 1.951
------------------------------------------------------------------------------
. tobit lwf age married children education, ll(0) nolog
Tobit regression Number of obs = 2,000
Uncensored = 1,343
Limits: Lower = 0 Left-censored = 657
Upper = +inf Right-censored = 0
LR chi2(4) = 461.85
Prob > chi2 = 0.0000
Log likelihood = -3349.9685 Pseudo R2 = 0.0645
------------------------------------------------------------------------------
lwf | Coefficient Std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
age | 0.052 0.006 9.08 0.000 0.041 0.063
married | 0.484 0.104 4.68 0.000 0.281 0.687
children | 0.486 0.032 15.33 0.000 0.424 0.548
education | 0.115 0.015 7.62 0.000 0.085 0.145
_cons | -2.808 0.263 -10.67 0.000 -3.324 -2.291
-------------+----------------------------------------------------------------
var(e.lwf)| 3.507 0.150 3.225 3.814
------------------------------------------------------------------------------MLE 的估计结果与 Tobit 模型的估计结果完全相同。
4. 结束语
本文主要介绍了最大似然估计和条件最大似然估计的理论过程，并使用网格搜索对 MLE 中的数值计算进行了简单的示例，最后使用五个范例说明了在 Stata 中如何使用 ml 命令进行最大似然估计。
本文的范例主要涉及线性和非线性的单方程模型。实际上多方程模型同样可以使用 MLE，具体理论包括似无相关回归、完全信息极大似然估计和有限信息极大似然估计等。此外，ml 命令还可以进行更加复杂的设定，我们将在之后的推文中介绍。
5. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh mle probit tobit, m
安装最新版lianxh命令：
ssc install lianxh, replace
专题：Stata命令 Stata新命令：面板-LogitFE-ProbitFE 专题：Stata程序 Stata：数值求解极大值及 MLE 示例 极大似然估计 (MLE) 及 Stata 实现 专题：回归分析 多层级Tobit模型及Stata应用 专题：内生性-因果推断 Stata: 一文读懂 Tobit 模型 专题：交乘项-调节 Logit-Probit中的交乘项及边际效应图示 专题：Probit-Logit Logit-Probit：非线性模型中交互项的边际效应解读 详解 Logit/Probit 模型中的 completely determined 问题 二元选择模型：Probit 还是 Logit？ Stata：二元Probit模型 动态 Probit 模型及 Stata 实现
课程推荐：因果推断实用计量方法
主讲老师：邱嘉平教授
🍓 课程主页：https://gitee.com/lianxh/YGqjp
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。